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Abstract 

The parallel dynamics of the fully connected Blume-Emery-Griffiths neural network 
model is studied at zero temperature using a probabilistic approach. A recursive 
scheme is found determining the complete time evolution of the order parameters, 
taking into account all feedback correlations. It is based upon the evolution of the 
distribution of the local field, the structure of which is determined in detail. As an 
illustrative example explicit analytic formula are given for the first few time steps 
of the dynamics. Furthermore, equilibrium fixed-point equations are derived and 
compared with the thermodynamic approach. The analytic results find excellent 
confirmation in extensive numerical simulations. 

PACS : 87.10.+e, 02.50.+S, 64.60.Cn 
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1 Introduction 



Recently, an optimal Hamiltonian has been derived in the statistical mechanics approach to Q- 
state neural networks starting from the concept of mutual information [1],[2]. Optimal means 
that the best retrieval properties are guaranteed including, e.g., the largest retrieval overlap, 
loading capacity, basin of attraction, convergence time. For Q = 3 this Hamiltonian resembles 
the classical Blume-Emery-Griffiths (BEG) Hamiltonian in the sense that it contains both a 
bilinear and biquadratic term in the spins [3]. For a fully connected architecture it has been 
shown, using a thermodynamic replica approach that the maximal loading capacity for the 
BEG network is indeed bigger than the one for other three-state networks existing in the 
literature (Ising, Potts...) [2]. 
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of an asymmetric extremely diluted architecture where one knows that there are no feedback 
loops complicating the time evolution, this dynamics has been solved in closed form [1] showing 
a better retrieval quality than other diluted models for certain parameters of the system. For 
symmetric architectures - fully connected but also extremely diluted - however, the situation 
is much more complicated. From former experience with, e.g., QTsing models (see [4] and 
references therein) one knows that the dynamics is very non-trivial due to the feedback in 
the system [5]. This feedback causes the appearance of discrete noise, besides Gaussian noise, 
involving the neurons at all previous time steps and prevents a closed-form solution. 

In this work we generalize the probabilistic approach that has been developed for the Hopfield 
model [6] , [7] and Q-Ising model [8] in order to solve the dynamics for the fully connected BEG 
model at zero temperature. Thereby, we start from the time evolution of the distribution of 
the local field, instead of working directly with the order parameters. We study the structure 
of this distribution in detail and, using this knowledge, we develop a recursive scheme in order 
to calculate the relevant order parameters of the system, i.e., the main overlap, the neural 
activity, the activity overlap and the variance of the residual overlap at any time step. As an 
illustration we write out these expressions in detail for the first few time steps of the dynamics. 

Furthermore, by requiring the local field to be time-independent, implying that some cor- 
relations between its Gaussian and discrete noise parts are neglected, we derive fixed-point 
equations for the order parameters. They coincide with those derived via thermodynamical 
methods [2]. 

Finally we perform numerical simulations of a BEG network with N = 6000 neurons. They 
confirm the analytical results we have derived. 

The rest of this paper is organized as follows. In Section 2 we introduce the model, its dynamics 
and the relevant order parameters. In Section 3 we use the probabilistic approach in order 
to derive a recursive scheme for the evolution of the distribution of the local field, leading 
to recursion relations for the order parameters. Using this scheme, we explicitly calculate in 
Appendix A the order parameters and in Appendix B the local field for the first few time steps 
of the dynamics. In Section 4 we show the existence of a Lyapunov function at zero temperature 
and we discuss the evolution of the system to fixed-point attractors. Section 5 details the 
structure of the local field distribution, especially the appearance of gaps. The analytic results 
are compared with numerical simulations in Section 6. Some concluding remarks are given in 
Section 7. 



2 The model 



Consider a neural network consisting of N neurons which can take values a^i = 1, . . . ,N 
from the discrete set S = { — 1,0, +1}. The p = aN patterns to be stored in this network 
are supposed to be a collection of independent and identically distributed random variables 
(i.i.d.r.v.), /i — 1, . . . ,p with a probability distribution 

p(tf) = \m - 1) + f w + 1) + (i - a)m) (i) 
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Given the network configuration at time t, cr N (t) = {crj(t)}, j = 1, ... ,JV, the following 
dynamics is considered. The configuration <tjv(0) is chosen as input. At zero temperature all 
neurons are updated in parallel according to the rule 

<Ji(t) -> (Ti(t + 1) = s' : minej[s|o-jv(0] = e i[ s Viv(0] ( 3 ) 

with s' G 5. We remark that this rule is the zero temperature limit of the stochastic parallel 
spin-flip dynamics defined by the transition probabilities 

Pr( gi (t + 1) = s' 6 S\* N (t)) = ^ V/fL ■ < 4 ' 

Here the energy potential ej[s|<Tjv(i)] is defined by 

€i(s\cr N (t)) = -shi(o- N (t)) - s 2 9i(cr N (t)) , (5) 
where the following local fields in neuron % carry all the information 

h N ,i(t) = E -W*), M*) = E (6) 

with the obvious shorthand notation for the local fields. The synaptic couplings Jij and 
are of the Hebb-type 

J* = 4* E tftf, = ^ E « (7) 



with 



These are the local fields entering in the BEG model [1]. The updating rule (3) is equivalent 
to using a gain function 

<7i(t + l) =g(h Nti (t),e N ,i(t)) = siga{h N>i {t))Q{\h N<i {t)\ +6 N>i (t)) (9) 

with 6 the Heaviside function. 

The order parameters of this system have been obtained starting form the mutual information 
as a measure for the retrieval quality of the system [1], [2]. They are the retrieval overlap, the 
activity overlap, and the neural activity 

™&(*) = ^E#*(*), <W = ^E(en 2 (^W) 2 , ?^)4W)) 2 - ( 10 ) 
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activity overlap 



3 Recursive dynamical scheme 

In networks with symmetric couplings it is known that non-trivial correlations occur, even at 
zero temperature, which become increasingly tedious to evaluate [5]. 

On the basis of the probabilistic approach (see, e.g., [6], [7], [8]) used successfully before ([4] 
and references therein) we develop in this section a recursive dynamical scheme in order to 
study the time evolution of the distribution of the local fields hi(t) and $i(t). This allows us 
to write down recursion relations determining the full time evolution of the order parameters 
(10)-(11) of the BEG network model. 

Suppose that the initial configuration of the network {o"j(0)} is a collection of i.i.d.r.v. with 
mean E[<7j(0)] = and variance Var[<7j(0)] = go an d correlated with only one pattern which 
we choose, without loss of generality, to be the first one 

E[£V,-(0)] = S^mla, m\ > 0, E[t#<7?(0)] = S^l* . (12) 

By the law of large numbers (LLN) eqs. (10)-(11) and (12) determine the order parameters 
m]y(0), gjv(0) and 1^(0) at t = in the limit iV — > oo. 

Next, we want to apply standard signal-to-noise techniques (see, e.g, [6], [8]) to both the local 
fields fojv,i(0) and 6*^(0) at t — 0. Starting from their definitions we find 

£l ? X(0) + Af(0,^) (13) 

em = jim Uk(p) - Uvi?am + jr EEl?»)) 

£, V(0)+A ,(o,-^), (14 ) 

where the convergence is in distribution [9]. The quantity Af(0,d) represents a Gaussian ran- 
dom variable with mean and variance d. 

The key question is then how these quantities evolve in time under the parallel dynamics 
specified before. For a general time step we find from eq. (9) and the LLN in the limit N — > oo 
for the order parameters (eqs.(lO)-(ll)) 



4 



m^ + l^-^MO,^)))) (15) 
^ + 1) = ((CV(^W,^W))) (16) 

i^t + i^^lhiWMt)))), (17) 

where hi(t) = linijv^oc h N>i {t) (with an analogous formula for 9i(t)), and where the convergence 
is in probability. In the above ((•)) denotes the average both over the distribution of the {£f } 
(and hence {^f}) and the {<7j(0)}. Note that the average over the latter is hidden in an 
average over the local field through the updating rule (9). From the work on symmetric Q- 
Ising networks [4], [10] we know that due to the correlations we have to study carefully the 
influence of non-condensed (/i > 1) patterns in the time evolution of the system, expressed by 
the variance of the residual overlaps, in our case in both the local fields. The latter are defined 
as 



r"(f) = Mm r%(t) = lim -±= £^(t), p > 1 (18) 

N->oo N^oo a 2y/]\[ ^ 

^(^lim s%(t) = hm -L^^®' H>1 (19) 

N^oo N—*oo _/Y ~~ 

where the limit iV — > oo of r^(0) and s%(0) is given by the Gaussian random variable in eqs. 
(13) and (14). At this point we want to remark that the choice of the initial configurations 
assures the independence of r M (0) and s M (0) as can be seen by calculating the characteristic 
function E[exp(ixr^(0) + iys%(0))]. 

The further aim of this section is then to calculate the distribution of the local fields and the 
order parameters as a function of time. 

We start by rewriting the local fields (6) at time t in the following way 

hN.it) = h}m\(t) + J- £ £ £f^(0 - % t (t) 
a a iV n>i j a 

= klm^t) - % t (t) + J= £ £W) ( 20 ) 
a a ViV M>1 

1 n 

= tf&(*) - 77T^-T^(^) + E «(0 • (21) 
a(l-a) V-/V M >i 

From a technical point of view the explicit addition and subtraction of the (Ti(t) (af(t)) term is 
convenient in order to treat all indices in the sum over j on equal footing, which is important 
to take into account all possible feedback loops. 

At this point several remarks are in order. Since the neuronal states {<Jj(t)}, for t > 0, are 
not i.i.d.r.v., the central limit theorem (CLT) can not be applied directly to the residual 
overlap r^(t) and s%(t). Furthermore, the set of aN variables an( l {^ s at(0}^ 

are not independent because the r v N {t) respectively s v N {t),u ^ ji are weakly dependent on 
the £f respectively rjf. Indeed, after applying the dynamics, the <Ji(t) (of(t)) and the £f (r?f ) 
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dependence gives rise to a macroscopic contribution after summing and taking the limit iV — > 
oo. Therefore, we follow a procedure similar to the one used in the Q-Ising model [4], [10] by 
isolating in the local fields precisely the contributions arising from these dependences. 

In order to do so we rewrite the residual overlaps as 



A" ' 



TV 
=? 



(22) 
(23) 



with obvious notation. In these expressions we have extracted the contribution of the \x term 
out of the local fields such that the modified local fields h%^{t) and Q%i(t) are only weakly 
dependent on £f and respectively, whereas h^^t) and O^^t) depend strongly on them. 

Next, we want to find the most important terms in (22) and (23) in the limit N — > oo. 
Therefore, we consider the characteristic function E[exp(ixr^(t + 1) + iys^{t + 1))] using (22) 
and (23), up to order 0(N~ 3 / 2 ). We then expand the gain function around the modified local 
fields. After some calculation we obtain in the limit A" — > oo 



lim E[exp(ixr^(t + 1) + iys%(t + 1))] 

N^oo 



exp 



ix Xh (t)r»(t) - * 2 ^±^ + iyxe(t)s»(t) - q{t + 1} 



2a(l - a) 



(24) 



with Xh(t) and Xe{t) the "susceptibilities" corresponding to the fields hi(t) and 9i(t) and given 
by 



Mt) = (( 1 1» = ((; /I * C <*««>< w> (IIJ » 

= | ((jT <*>«,(,) (<%«)(<?))) + (1 " <>)»(*) (25) 

„)) = 5(1^0 ((yl^^^ + «»(-*»>) ' (26) 

In these expressions, p^Jh) and Pqu\(0) are the probability densities of the modified local 

fields, /i and 0. They are the integrations of the joint distribution p~ h ^ §^(h, 6) over the h and 
9 values, e.g. p%^{h) = J dOp^ §^j(h, 6). (See Section 5 for more details). From the expansion 
(24) we see that the local fields h N:i (t) and 9 Nji (t) are independent up to the order 0(N~ 3 ^ 2 ) 
since r^(t) and s^(t) are as well. 

Identifying terms we then get 



Xe{t) 




r^(t + l) = ^(t)+ X h{t)r^t) (27) 
s' i (t + l) = s*(t)+xe(t)s*(t), (28) 



where 
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i(l — a) 

We remark that this calculation also shows us that r M (t) and s M (t) are independent for all 
times. 

In this way we obtain in the limit iV — > oo from eqs.(20) and (21) 



hi(t + 1) = -£rn\t + 1) + - hlm\t) + -a,(t)l 

a I a a J 

+A/-(o,^g(t+l)) 
fc(f + 1) = ^(f + 1) + X e(t){9i(t) - r,]l\t) + ^ 2 (t) 



(31) 



+a/Yo, 



V ' a 2 1 - a 



??(*+!)■ 



(32) 



From this it is clear that the local fields at time t + 1 consist out of a discrete part and a 
normally distributed part, viz. 



hi(t + 1) = Mi(* + 1) + A/"(o, V(t + 1) 
0<(* + 1) = L<(* + 1) + A/"(o, W(* + 1)) , 

where 



(33) 
(34) 



Mi(t + l)= X h(t) 
L i (t + l)= X e(t) 



M i (t)-^m 1 (t) + -a i (t) 
a a 



a 



Li(t) - v}l\t) + 



a 



a(l — a 



and 



V(t + l) = aaD(t + l), W(t + 1) 



a 



a(l — a) 



(35) 
(36) 

(37) 



with D(t + 1) and E(t + 1) the variances of the residual overlaps, r^(t + 1) and s^(t + 1), 
satisfying the recursion relations 



D(t + 1) = + xl(t)D(t) + 2 Xfc (t)Cov[P(t), r"(f)] 

E(t + 1) = + ^(0^(0 + 2 X ,(f)Cov[S"(f), *"(*)] • 

all — a 



(38) 
(39) 



We still have to determine p~ h( ^(h) and p§^(9) in (25) and (26). We know that the quantities 
Mi(t) and Li(t) consist out of a signal term and a discrete noise term, viz. 



7 



M i (t) = *-m 1 (t) + 'E- Uxh(s) atf) 



t'=0 

t-i 



L l {t) = V ]l\t) + Y J 



a 



, =0 a(l - a) 



t-i 



II Xe(s) 



s=t> 



(40) 
(41) 



The evolution equation tells that the <Jj(t') and uf{t') can be written in terms of the hi(t' — 1) 
and 9i(t' — 1) such that the second terms in the expressions above are the sum of correlated 
variables. Furthermore, these are also correlated through the dynamics with the normally 
distributed part of the local fields. So the local fields can be considered as a transformation 
of a set of correlated variables x = {x s }, y = {y s }, s = 1,2,... , t — 2, t which we choose to 
normalise. Then we arrive at the following expression for the probability densities of the local 
fields 



J im Ph»(t)( h ) = Phi(t)(h) = / ( TT dx s dy s ) dxtdyt—^ —= 

N ^°° A ) J \=o 1 ^&et(2nC h ) ^/det(27rC e ) 

exp ( - ixC h - x x r - \yC^y T )5(h - M t (t) - y/vftxt) 

r l ~ 2 1 
lim p#. (t) (0) = p 9 . {t) (0)= / (]Jdx s dy s )dx t dy t = 

1 J s=o \/det(^ 



(42) 



det(27rC /l ) ^det(2irC g ) 

exp ( - \*C?iF - \yC^y T )5(e - L t (t) - Jw(t)y t ) , 

where the correlations matrices Cy L and Cg are given by 



(43) 



(44) 



Together with eqs. (15)-(17) the equations (27)-(28),(35)-(39),(42)-(44) form an exact recursive 
scheme in order to obtain the order parameters of the system. The practical difficulty that 
remains is the explicit correlations in the network at different time steps. As an illustration 
we calculate the first three time steps in Appendix A. 



4 Fixed-point equations 



A second type of results can be obtained by requiring through the recursion relations (35)- (39) 
that the local fields become time- independent. This means that most of the discrete noise part 
is neglected. We show that this procedure leads to the same fixed-point equations as those 
recently found from a replica symmetric thermodynamic approach [2]. 

First, for the BEG-model one can show that 

Hit) = -£ [h^m + ei(t)5?(t) , (45) 
i=i v 7 
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that 

ei [^(t)|o-(t)] = mine i [s| < r(t)] (46) 

is a Lyapunov function for zero temperature. The proof is completely analogous to the argu- 
mentation used in [11] and [12]. The choice of {&i(t)} implies through the updating rule (3) that 
<Ji{t + 1) = o~i(t). For finite N, H(t) is bounded from below implying that H(t + 1) — H(t) = 
after finitely many time steps. This can be realized for Gi{t + 2) = <Ji(t) for all % and, hence, 
both two-cycles and fixed-points satisfy this condition. We only study fixed-points. 

Next, we start by eliminating the time dependence in the evolution equations for the local 
fields (31)-(32). This leads to 



1 1 . r( ol \ a , 

h l = -i l m + - M[ 0, —q + -r) h ai 47 

a 1 — Xh \ a J a 

| / Oi \ Oi 

0i = Vil + M 0, — -q + — , m a] , (48) 

1 — xe V a 2 (l — a) 2 J a(l — a) 

where from now on we forget about the pattern index 1 and where we have defined 

rix = T^—, x = h,9. (49) 
1 - Xx 

This means that out of the discrete part of the local field distributions, i.e., M^t) (Lj(t)), only 
the <7j(i — 1) (of (t — 1)) term is kept besides, of course, the signal terms. These expressions 

consist out of two parts: A normally distributed part, hi = Afi^im/a,aq/a 2 (l — Xh) 2 ) and 



the analogous formula for $i, and some discrete noise part. Employing these expressions in the 
updating rule one finds 

/- a ~ a 2 \ 

(Ti = g\hi + -r) h Gi, Oi + — -rio^ . (50) 

V a a(l — a) J 

This is a self-consistent equation in <jj which, in general, admits more than one solution. This 
type of equation has been solved in the case of analog neural networks with continuous time 
dynamics [13] and in the case of Q-Ising neural networks [4], [10] using a Maxwell construction. 
Here we follow the same line of reasoning for the joint probability distribution of the local 
fields in the (h, 0)-plane (see fig. 1) leading to a unique solution 

ai = g{h h e)j =sign(^)©(N +K + a) , (51) 

where 

= Ta'» + 2^0^) *' < 52 > 

Using the definition of the order parameters (see (10), (11)) in the limit iV — > oo one finds in 
the fixed point, dropping the index i 
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From (27)-(28), (38)-(39) and (25)-(26) it is clear that 

n _ g p - g 

a 3 (l-X,) 2 ' a(l-a)(l- X .) 2 

with 



= y = 1 = = (( f Dz [ Dy y g 2 (-£m + - z, rjl + — ry) )) . 

yja(l - a)aE W ^ V « a(l ~ Xfc) a(l - a)(l - xe) 7 // 

(58) 

These equations are the same as the fixed-point equations derived from a replica-symmetric 
mean-field theory treatment [2]. 
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Fig. 1. The Maxwell construction. To the right of the short-dashed line the solution a = 1 exists, to 
the left of the long-dashed line, a = —1 is a solution, and under the dotted line, a = exists. The 
thick full line shows the unique solution, dividing the local fields space in three parts. The point P is 
given by P = (o, -A) . 




5 Local field distribution 



It is interesting to study the distribution of the local fields, the main ingredients in our dy- 
namical scheme. First, we look at the stationary distribution. Since the Maxwell construction 
we used is discussed in the plane (h,9), we want to find the joint-distribution for the local 
fields, Poo(M) = ph(oo),9(oo)(h,9) defined by 
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Pco (h, 9) = ^J DzDyS [h - -Jm - - Vh a - —J^) 

xs(o -rjl- -7^—,Vev 2 ~ t^—J 2n a<1 . 2 ~yW) , (59) 
V a(l - a) 1 - xe y a (1 _ a ) > 

where $(cr) is obtained from the updating rule after the Maxwell construction (see eq. (51)) 



$(tr) = e(- -0 



cr,0 



2a 2a(l-a) 
This leads to 

Poo(K 9) = P+1 (h, 0)Q(h - % h )Q{h + 9 - A) + p^h, 9)Q(-h - % h )Q{-h + 9 - A) 

+ A) (/ l ,0)e(-|/ l |-0-A), (61) 

where 



(l- Xh )(l- Xe )(l- a y 
2nqa 



exp i - -- 



1 (/i 



aq 



x exp i - -- 



2\2. 



i (* - - ^ W) 



aq 



(1-Xe) 2 a2(l-a)2 



(62) 



Analyzing these expressions we see that the distribution Poo(h, 9) shows a gap. In fig. 2 we show 
this gap structure which depends, of course, on the specific values of the physical parameters 
a, a, Xh, Xe of the system. The important points bordering these gaps are given by P + i = 
O,A-»,P_ 1 = (->,A-»,P = (0,-A). 





G 


e 


p 

+1 


Py 

A? 


e 




h 



Fig. 2. The gap structure of poo(h,8). The coordinates of the points P+i,P_i,Po are given in the 
text. The integral p-ti in (61) is only different from zero in the region to the right (left) of the line 
on which P±\ lies; po exists only below the line on which Pq lies. 
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projections on the 9 (or h) axis 



pt=oo{h\0 ) = 

(d\h ) = 



Pt 



Poo (Mo) 

S?ooPoo(h,e )dh 

Poo(h ,9) 

JZoPoc(h , 



(63) 
(64) 



Finally, starting from (42)-(43) one can write down expressions for Ph{t)(h), Pe(t)(Q) f° r the 
first time steps and calculate both the joint probability and its projections from it. This is 
illustrated in Appendix B. All these projections will be compared with numerical simulations 
in the next Section. 



6 Numerical results and simulations 



The equations derived in Sections 3-5 and Appendices A and B have been studied numerically 
and have been compared with simulations for systems up to N = 6000 neurons averaged over 
500 runs. 

We start with the remark that the initial conditions are not independent because positivity of 
the relevant probabilities implies 



g > am , 



> Iq, and for a > go : 



for a < g : Zq < 



% 
a 

1 - go 
1 — a 



(65) 



The phase diagram of the fully connected BEG neural network has been discussed in [2] using 
a replica-symmetric mean-field theory. From that work we see that for uniformly distributed 
patterns the critical capacity at zero temperature is 0.091. 

The first point we would like to examine is whether the recursive dynamical scheme we have 
derived is confirmed by simulations. This is illustrated by some typical results in fig. 3 showing 
the order parameters as a function of a for uniform patterns and m = 0.6, l = 0.6, g = 0.5. 
We see that the theoretical results, given by the explicit formula in Appendix A, and the 
simulations agree very well over the whole range of a's. Furthermore, we learn that in the 
retrieval regime the first time steps of the dynamics give us already a reasonable estimate for 
the critical capacity especially through the order parameter /. This is also the case for the 
other values of a, a = 0.02, 0.05, 0.08, we have considered. 

These findings are confirmed by some typical (m(t), l(t)) flow diagrams. In figs. 4 and 5 we show 
the results for uniform patterns and two values of a in the retrieval region, a = 0.015, 0.08, 
giving us a good idea of the basin of attraction. Remark that to the right of the dotted line 
we cannot start initially because of the condition (65). At later times we can enter this region 
because q(t) changes from its initial value g(0) = a = 2/3. The dashed line in the figures 4a 
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q(t) 




0.1 

a 



Fig. 3. Order parameters m(t), l(t) and q(t) as a function of the capacity a for the first three time 
steps and a = 2/3, mo = 0.6, = 0.6, qo = 0.5. Theoretical results (solid lines) versus simulations 
with N = 6000 (time 1, 2 and 3 given by the plus symbol, times symbol respectively circles) are 
shown. 

and 5a indicates the border of the basin of attraction. To have an idea about the accuracy 
of this basin boundary we also show in figs 4b and 5b the percentage of runs going to the 
attractor on the line mo = Iq, starting from the (0,0)-point. The basin boundary is drawn 
joining the starting points of the flow lines reaching the attractor with a percentage lying 
between 45% and 55%, as visualized by the two parallel dashed lines. As expected, the basin 
of attraction shrinks for increasing a, due to the appearance of other thermodynamically stable 
states (spin-glass states), and the error for defining the basin boundary becomes bigger. For 
comparable values of the relevant system parameters for the Q-Ising model one can verify (see 
[10]) that the basin of attraction of the latter is smaller. 
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0.2 0.4 0.6 0.8 1 0.04 0.08 0.12 

m(t) m(0)=l(0) 



Fig. 4. Flow diagram for uniform patterns a = 2/3 and capacity a = 0.015. In figure 4a, the 
long-dashed line shows the basin boundary. The parallel dashed lines in figure 4b indicate the error 
bounds. 



(a) (b) 




0.2 0.4 0.6 0.8 1 0^2 0.4 0.6 0.8 

m(t) m(0)=l(0) 
Fig. 5. As in fig. 4 for a = 0.08. 



The next point we look at is the appearance of a so-called quadrupolar state (m = 0, 1 ^ 0). 
Extended simulations did not show any such stable state at zero temperature in the fully 
connected model. This is different from the findings in the diluted model where a quadrupolar 
state is predicted [1]. However, recent discussions indicate that also for this diluted model the 
quadrupolar state is not stable at temperature zero [14]. 

Finally, we consider the distributions of the local field, an important ingredient in our dynami- 
cal scheme. We have investigated them numerically using the fixed-point equations mentioned 
before and compared them with numerical simulations. Some typical results are shown in 
figs. 6 for uniform patterns. In fig. 6a we show the joint distribution of the local fields, for 
£ = 0, projected on the /i-axis, fixing 9 = —0.6, in the retrieval region, a = 0.064, for 
m o — h — 0.5, q = 0.5, while fig. 6b represents this distribution in the non-retrieval region, 
a = 0.13, for m Q — l — 0.2, q = 0.5. The first time steps, given by the explicit formula in 
Appendix B, are in complete agreement with the numerical simulations. 

Concerning the gap structure we see that for the retrieval state there are, typically, small gaps 
in the equilibrium distribution. For small a the gaps are very narrow (see insets of fig. 6a). In 
the simulations these gaps show up very quickly. In order not to overload the figures we have 
plotted one intermediate result for the time step t = 100. For the non-retrieval state the gaps 



14 



cue ty piuaiiy iiiul.ii uiggci . -rvgcun 111 niiiiuidijiuiio tiic gajjo buuw ujj idunci l[u.il..i^ij<. 111 ng. i wc 

plot a 3-d picture of the equilibrium distribution for £ = in the spin-glass region, a = 0.013, 
for uniform patterns. The gaps are clearly visible. We recall that the theoretical equilibrium 
results coincide with the thermodynamic replica-symmetric solution and that it is expected 
that the gaps are reduced to one point for the exact solution. It is extremely difficult to find 
points touching the axis in the simulations because of the final size effects. Analogous results 
have been found for the Hopfield model, the Q-Ising model [15], [16], [17] and, first, in the 
infinite range spin- glass [18]. 



(a) (b) 




-1.5 -1 -0.5 0.5 1 1.5 -4 -3 -2 -1 1 2 3 4 



h h 

Fig. 6. Projections of the joint probability distribution of the local fields in the retrieval (a = 0.064, 
fig. 6a) and the spin-glass phase (a = 0.13, fig. 6b). Simulations for time steps (circles), 1 (squares) 
and 100 (plus symbol) are shown. For clarity we have not included time results in fig. 6b. Dotted, 
solid and dashed lines show the theoretical results for times 0, 1 and oo. 



7 Conclusions 

An evolution equation is derived for the distribution of the local field governing the parallel 
dynamics at zero temperature of BEG networks. All feedback correlations are taken into 
account. This distribution contains both a normally distributed part and a discrete part. 

Employing this evolution equation a general recursive scheme is developed allowing one to 
calculate the relevant order parameters of the system, i.e., the retrieval overlap, the activity 
overlap and the neural activity for any time step. This scheme has been worked out explicitly 
for the first three time steps of the dynamics. 

Under the condition that the local field becomes time- independent, meaning that some of the 
discrete noise is neglected, fixed-point equations are obtained for the order parameters. They 
agree with those obtained from a mean-field replica symmetric thermodynamic approach. The 
gap structure of the equilibrium local field distribution is examined. The gaps in the retrieval 
regime are much smaller than those in the non-retrieval regime. 

Extensive numerical simulations are performed for a system of 6000 neurons. They confirm the 
results obtained from the dynamical scheme both for the local fields and the order parameters. 
Furthermore, they illustrate that the first few time steps do give a reasonable estimate of the 
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P(h,e) 




Fig. 7. A 3-d plot of the joint probability for the local field distribution for £ = in the spin glass 
phase, a = 0.013, and uniform patterns, a = 2/3, at t — ► oo. 

critical capacity, especially through the activity overlap order parameter. Finally, flow diagrams 
indicate the size of the basin of attraction of the retrieval state as a function of the loading. 
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Appendix A 



Following the general recursive scheme developed in Section 3, evolution equations are derived 
for the first three time steps of the BEG fully connected network, taking into account all 
correlations. Our starting point is the set of equations for the order parameters (15)-(17) with 
the following initial conditions 

m 1 (0) = m , l\0)=l o , q(0)=q (66) 

HO) = klm +M(0, ), 0,(0) = r£lo+W, 2( T v? ) ( 67 ) 
a cr er(l — a) A 
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First time step 

We immediately get from (15)-(17) 

m(l) = ±((t J DzJ Dyg(K(*)Mv)) 
q(l) = t^J DzJ Dy g 2 (h' (z) , 9' (y] 
l{\) = ^r,jDzjDy g 2 (h' (z) , 6' (y) 

where Dx = dxexp(— x 2 j2)j\phx denotes the Gaussian measure and 

ti (z) = km + y/wjz, 9' (y) = rjl + Jw(0) y . 



(69) 
(70) 
(71) 



(72) 



From the definition of h^^t) and O^^t), we know that, when N — > oo, the elements in the 
pairs {ef,^(^(0),C(0))}i {£>,(0)}, K , <?W (0), 0f(O))} and K , ^(0)} are uncorrelated 
for /i ^ 1. Therefore, using the recursion relations (38) (39) we get 



£>(1) = + xl(0)D(0) + 2 Xfc (0)i2(l, 0) 

E(l) = -77^ + xl(0)E(0) + 2 X ,(0)S(1, 0) , 
a(l — a) 

where, in general, the correlation parameters are defined as 



(73) 
(74) 



R(t, = -^E^ - 1), 0(f - l))g{W - 1), 0(f - 1)) 



R(t,0) = -E <j(0)g(~h(t-l),9(t-l)) 



S(t,t') 
S(t,0) 



1 



a(l — a) 
1 



g^t-l),9(t-l))9 2 Ch(t'-l),0(t'-l)) 



a 2 (0)g 2 (h(t-l),9(t-l)) 



a(l — a) 

leading, in the limit N — > oo, to the following formula for the first time step 



(75) 
(76) 
(77) 

(78) 



R(1,0) = ±((*(0) J DzJ Dy g(h' (z),9' i 
S(l,0) = -^—- ) ((^a 2 (0) jDzjDy g 2 (h' (z),9' (y))^ . 

Since at time zero there are no correlations yet 



(79) 
(80) 
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X„(0) = — ^ {{J DzJ Dyz g(h' (z), ff (»))ft (81) 

*« (0) = JTzhmSi ((/ Dz J Dyy 9 ^ ^ » ' (82) 



a(l-a)yJW(0) 
Second time step 



First, we need the distribution of the local fiels at time t — 1. This follows immediately from 
(35)-(36) and (37) 

fci(l) = ^m(l) + - Xfe (0)^(0) + A/"(0, V(l)) (83) 

e i (l)=r ]l l(l) + ——-Xo(0)aKO)+mO,W(l)) . (84) 

(X 1 J. ft ) 

These results allow us to write down the order parameters at time step 2: 

m(2) = i^ J DzJ Dy g(K(*)Mv)))) (85) 

?(2) = ((/l>* / £>y <7 2 (/*i(*),^(l/)))) (86) 

Z(2) = ^ /^/ Dy ffcWMvj))) , (87) 



where 



= -$m(l) + -X,(0)^(0) + y/vO)z (88) 



a 



e[(y)= V l(l) + __- Xff (0)(7?(0) + (89) 
a(l — a) v 



and 



(90) 



M(1) = ad-aljm (d DZ I Dy ' ■ (51) 

The calculation of the variance of the residual overlap needs some more work. From the 
recursion relations (38)-(39) one finds 

D(2) = ^ + xl(l)D(l) + 2 Xh (l) (R(2, 1) + Xh (0)R(2, 0)) (92) 

E{2) = + xlmil) + 2X8(1) (5(2, 1) + **(0)S(2, 0)) , (93) 
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i2(2, 0) = 1 ^r(O) £>y fl (/ii(«),^(y)) 

^(2,0) = ^ T ^((,x 2 (0) g 2 (K(z),e[(y)) 



(94) 
(95) 



To obtain R{2, 1) and 5(2, 1) we remark that the local fields at time steps and 1 are corre- 
lated. The correlation coefficients of their normally distributed part, viz. 



E[(fo(t) - M(t))(h(f) - M(t'))} 
Ph{t, t ) = ■==—== 

E[(9(t) - L(t))(h9(t') - L(t>))} 
Pe{t,t) = 



'W(t)y/W(V) 

is found using the recursion formula (83)- (84) 

W, ) + D(Q)xM 



'D(Q)D(1) 

Employing all this in eqs. (75) and (77) we arrive at 



S(l,0 ) + E(0)xe(0) 
E(0)E(1) 



i2(2,l) = l^/ D^\z,s) J D^°(y,t) g (h[(z) , e[(y)) g(h' (s),6> 



(96) 
(97) 



(98) 



(99) 



5(2,1) 



a(l — a) 

where the joint distribution Du^ ,b (z,y) equals 

dzdy 



Dcoj : (z,s) J Du\% y t) g 2 (h[(z) , O'^y)) g 2 (h' (s), ^(f)))) ,(100) 



DoJ a /(z,y) = 



2n^l -p? x (a,b) 



cxp 



z 2 -2zyp x (a, b) + y 2 
2(1 -pl(a,b)) 



(101) 



Third time step 



We start by writing down the distribution of the local fiels at t — 2 



h l (2) = ^m(2) + - Xh (2) 

(X (X 



^(1) + Xh(0)<7 i (0) 



+ M(0,V(2)) 



e l (2)= Vl i(2) + 



a 



a(l — a) 



a 2 (l) + X e(0)a 2 



+ Af(0,W(2)) . 



(102) 
(103) 



In order to write down the expressions for the order parameters starting from (15)-(17) the 
average has to be taken over the Gaussian noise, o"j(0) and <7j(l). The average over <7j(0) 
causes no difficulties because this initial configuration is chosen randomly. The average over 
the Gaussian random noise variable appearing in hi(2), 0j(2), and <7j(l) is more tricky because, 
e.g., hi(2) and <7j(l) are correlated by the dynamics. However, the evolution equation tells us 
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instead of <7j(l). From the recursion relation (31)-(32) one finds for the relevant correlation 
coefficients 



Ph(2,0) 



P*(2,0) 



fl(2, 0) + g(l, O)xh(l) + g^(%(0! 

}/D(0)D(2) 
S(2, 0) + 5(1, 0) Xe (l) + E(0) X e(l)xe(0) 



(104) 



(105) 



y/E(0)E(2) 



Using this we get 



m(3) = i^ / £>"h ,0 (*,s) / Dul%,t)g {^,8,1)^,8,1))^ 
?(3) = ((/ ^°(z, S ) | L>o; e 2 ' (y,^ 2 (^(z, S ,0,^(y,s,t)))) 
J(3) = ^ij | £fc# (z,*) y Di4%,t)grM(z,8,t)My,8,t)ty 



(107) 



(106) 



(108) 



with the joint distributions as defined before (see, (101)) and 



1 Q/ / 

/i' 2 (*, s, t) = -im{2) + - Xh (l)[g(ti (s), 6' Q (t)) + Xh(0)<r(0)] + y/V(2) z 



(109) 



^(y, S , =17/(2) + — ^— x.(l)[/(^o(s), ^(0) + X*(0V 2 (0)] + v^(2) 1/ • (HO) 
a i J. d ) 



In the same way further time steps can be calculated at the price of more complicated algebraic 
expressions. 



Appendix B 

We calculate explicitly the projected joint distributions for the local fields for the first time 
steps. Starting from (42)-(43) we obtain 




(112) 



(HI) 



and 
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Ph(i)(h) 



1 

^2nW(l) 



exp 



exp 



2V{1) ~ J 

-(fl-^a)-^^(oy 2 (o)) 2 

2W(1) 



}■ 



(113) 
(114) 



To construct the joint probability for t = we just have to take the product of (111) and 
(112). In order to find the joint probability for t — 1 we have to average the product of (113) 
and (114)over cr(0). This leads to 



Pt=i(h,9) = 



2nJV(l)W(l) 



x 



■ exp 



jh - Wi) - fx,(Q)) 2 -( g -^( 1 )-^yx e (o)) i 



2V(1) 



2W(1) 



^go J -(ft - H 1 ) + ^ ^(°)) 2 - - 

+ 2 exp 



+(1 - g ) exp 



2V(1) 
jh-jm{l)f 
2V{1) 



2W(l) 

2W(1) 



(115) 



The projected distributions are then obtained analogously to eqs. (63) and (64). For t = we 
find back (111) and (112) confirming again that the local fields at time zero are independent 
because no correlations are present yet. For the first time step we obtain 



Pt=i(h\9 ) 



27ry(l) 



go 

^exp 



2W(1) 



x 



+ 



X 



(exp j- 
(1 -g ) exp | 

{ 



■ (fe-£m(l)-f Xfc (0)) s 



2V(1) 
-(h-im(l)) 2 



+ exp 



2V(1) 



exp 



go exp 



2V(1) 

-(Q - V l(l) - ^jXeW) 
2W(1) 



2W(1) 



+ (1 - g ) exp 



-(6>p-^(l)) s 
2W(1) 



(116) 



and finally 
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Pt=i(0\h o ) 



2ttW(1) 



exp 



2W(1) 



x exp 



-(h Q - jm(l) - *Xh(0))'' 
2V(1) 



exp 



-(ftp - |m(l) + f x fe (0)) S 



2F(1) 



+ (1 - g ) exp \ ;t77777 r exp ' 



x 



<?0 

yexp 



2V(1) 

_ (/i0 _ | m(l) _ 2^(0)) = 



+ (1 - g ) exp 



2V(1) 

-(/iO-Hl)) 2 --! 



2F(1) 



2W(1) 
yexp 



-(hp - |m(l) + g Xfc (0)) s 
2V(1) 



;ii7) 
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